Stability of localized wave fronts in bistable systems 
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Localized wave fronts are a fundamental feature of biological systems from cell biology to ecology. 
Here, we study a broad class of bistable models subject to self-activation, degradation and spatially 
inhomogeneous activating agents. We determine the conditions under which wave-front localization 
is possible and analyze the stability thereof with respect to extrinsic perturbations and internal noise. 
It is found that stability is enhanced upon regulating a positional signal and, surprisingly, also for a 
low degree of binding cooperativity. We further show a contrasting impact of self-activation to the 
stability of these two sources of destabilization. 
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Bistable systems are ubiquitous in nature. For exam- 
ple, genetic switches are bistable systems that store the 
activation state of a gene [TJ[^. In population dynamics, 
a minimum population size is often needed to establish 
a stable population [3] . The spatial versions of such sys- 
tems admit traveling wave solutions, e.g. describing the 
outbreak of viruses or the colonization of territory (4HZ] . 
If bistable systems are subject to external spatial gra- 
dients, traveling waves may localize in restricted spatial 
domains [BHlOj. An important example arises in early 
embryogenesis of Drosophila melanogaster where mater- 
nal morphogen gradients provide positional information 
for gene regulation [TTI - fH] . The morphogen Bicoid is 
present as a monotonically decreasing concentration in 
the embryo and controls the step-like activation of the 
gene Hunchback, which also enhances its own activation 
effectively producing a bistable system. The exact po- 
sition of the Hunchback front is pivotal to the embryo's 
fate [13]. Hence, the front's stability to extrinsic pertur- 
bations or internal noise is paramount. Wave localization 
and the stability of the front also play an important role 
in other contexts. In ecology, birth rates may have spatial 
dependence, e.g. due to spatial variance in temperature 
or resource availability |15l [TB] . The localized boundaries 
between species are subject to large fluctuations due to 
the low number of individuals in the boundary region. 
This may eventually lead to the extinction of one of the 
species due to demographic stochasticity |17| . Last, in 
bio-technological applications, this mechanism might be 
used to create localized fronts of proteins . 

Motivated by these processes, we investigate a broad 
class of bistable diffusion-reaction models with reaction 
terms comprising self-activation, external activation, and 
degradation. While self-activation and degradation are 
assumed to be spatially uniform, the external activation 
is taken to be position-dependent. We consider two qual- 
itatively different types of external gradients and deter- 
mine the parameter range for which wave localization is 
possible. Moreover, we ask how stable these localized 
fronts are with respect to extrinsic and intrinsic noise. 



and we determine optimality conditions minimizing the 
front's susceptibility to such perturbations. 

Specifically, we consider a one-dimensional system 
where diffusing particles are subject to three types 
of reactions: First, there are gain processes with a 
concentration-dependent rate that accounts for self- 
activation in gene regulatory systems or reproduction in 
population dynamics. Typically, these rates are small 
for low concentrations, then rise and finally saturate at 
high concentrations. In populations dynamics, this is 
referred to as the strong Allee effect [31[B]. In gene reg- 
ulation, it can be due to cooperative transcription factor 
binding to a gene promoter. A common choice for the 
overall reaction rate is krK^gia) with the Hill function 
i?"^^ (a) = ^„'^^„ , kr the maximum intrinsic production 
rate, and a the particle concentration. The Hill coef- 
ficient n measures the degree of cooperative binding in 
the promoter region, or, in ecology, the strength of an 
Allee effect. Second, we account for loss processes, where 
particles vanish with a certain rate A. Third, in addi- 
tion to self-activation, there may also be external sources 
for particle production. Here, we are interested in sys- 
tems where this source is position-dependent and char- 
acterized by the overall rate kMM{x). The prefactor kjyj 
denotes the maximum rate of external activation, and 
M{x) is a monotonically decreasing positive density pro- 
file with normalization M(0) = 1. In the simplest case, 
where the profile results from a source-degradation dy- 
namics HH HO], it is exponential M{x) = e'^/^ with the 
decay length ^, cf. Fig. [l](a). Prominent examples are 
the concentration profile of Bicoid in Drosophila [19 and 
temperature or nutrient gradients in population dynam- 
ics [21 . Since the production of Hunchback by Bicoid 
is mediated by cooperative binding, the profile M(x) en- 
tering the overall production rate is commonly described 
by M{x) - R^^{e-^/^) The exponentially decaying 
signal induced by maternal source-degradation dynamics 
serves as an input to the gene regulation system. The 
latter is described by a Hill coefficient m typically in the 
range from 1 to 5, and an activation threshold /g. 
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Figure 1. (Color online) (a) Two types of gradients: exponential decay (dashed line) and a sigmoidal profile ensuing from 
regulating an exponentially decaying input (solid line), (b) The potential for different values of the front position q. The 
sliding ball analogy states that the front localizes where the two maximum values of the potential are equal, (c) Sketch of 
the bifurcation diagram and traveling wave solution of Eq. ([TJ . Blue lines denote stable solutions while the dashed (red) line 
corresponds to the unstable branch. Wave fronts (black lines and shaded area) penetrating the bistable region slow down and 
eventually come to rest at a stable fixed point of the front dynamics, (d) Phase diagrams of possible parameter values allowing 
wave localization. The parameter range of wave localization increases with the Hill coefficient n. 



In the limit of a large system size, fluctuations are 
of minor importance and the spatio-temporal dynamics 
is then aptly described by a reaction-diffusion equation, 
which in dimensionless form reads 



dtu = f{u,x) +dxxU. 



(1) 



Here f{u,x) = rR^^^{u) + M{x) — u comprises self- 
activation, external activation and degradation. Con- 
centration w, time t, and space x are measured in units 
of kM/K l/'^ and ^J~DJX, respectively. The ratio r = 
krjku denotes the relative amplitude of self-activation 
and external activation mediated through M{x). 

Traveling wave solutions of Eq. ([T|) may be localized 
due to the combined effect of spatially varying external 
sources and bistability [5fllO|. The basic mechanisms can 
be best understood in terms of the well-known sliding 
ball analogy |23], which here is complicated by the fact 
that the reaction term is space-dependent. Since in most 
biological situations a steep profile in u is induced by a 
smooth external profile Mix), we may assume a separa- 
tion of length scales ^ ^ ^J~DJ\ and ^ much smaller than 
the system size. Then one can make a generalized trav- 
eling wave ansatz V = lJ(x — q{t),y), where x is a fast 
varying variable describing changes in the concentration 
profile, y = x/^ is a slowly varying variable describing 
changes in the external profile M{x), and q{t) denotes 
the front position. To leading order this gives 

^qdxU^dxxU + duV{U,y)+Oir^), (2) 

which may be interpreted as a force balance for a par- 
ticle (sliding ball) with mass 1, friction q and potential 
V{u,y) — J f{u,y)du. Importantly, the potential para- 
metrically depends on y, see Fig. [ijb). For parameter 
regimes where V has two maxima at u^{x) and u~{x), 
and a local minimum at u*(a;), the velocity q must be 
chosen such that the sliding ball starting from the up- 
per branch u"*" ends up at the lower branch . The 



front speed is proportional to the difference between the 
two maxima of V{u, y) and becomes zero if the condition 

l^V{y) = /(m, y)Au — Q is satisfied. More quantita- 
tively, following standard steps |23fl25 ]. one finds 



^V{q) 



j-oo\-^xU{x - q,y)f dx 



(3) 



where U is the traveling wave solution. The denomina- 
tor roughly equals the maximum steepness of the front 
profile, and implies that steep fronts move slower [23J. 

In our class of models, a single branch of stable so- 
lutions at high concentrations typically undergoes a fold 
bifurcation for growing x, where the system is bistable on 
a confined spatial interval, see Fig. [l](c). For large X, a 
single branch at low concentrations remains. Within the 
bistable regime, the velocity c{q) may change sign and 
thereby lead to a localization of the traveling wave front. 

We first determine the localization position go of the 
front from /^V{qQ) = 0. Approximations for u^{x) can 
be obtained by expanding / as Taylor or Laurant series: 
u'ix) = M{x) + 0{u") and u+{x) = M{x) + r + 0(u-"'). 
For a given external profile M{x), the potential reads 
V{u,x) = -u[u/2 - M{x)-r + rF{u"/u^)], where 
F{z) = 2^1 (1, l/n ,1 + 1/n , —z) and 2F1 signifies the 
hypergeometric function. Keeping the dominant terms 
of M{x) in AV we then obtain an expression for Mq = 
Af ((Jo) determining the localization position q^: 
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This is well approximated by a linear function of the form 
g{n) ■ {uq — r/2) and converges to uq — -^r for n ^ 00. For 
exponentially decaying gradients, the equilibrium front 
position is then given by go = C In Mq- For sigmoidal gra- 
dients, M{x) = k ■ i?™(e~^/^) with dimensionless thresh- 



3 




(b) 1.0 

0.8 
o 0-6 



0.2 
0.0 

2 3 4 5 6 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 50 ^ 150 200 250 

r r Mo qg 

Figure 2. (Color online) Stability ^a, normalized to the steepness of the external profile, for (a) exponential and (b) sigmoidal 
[m = 5] external profiles; the Hill coefBcient for self-activation is n = 5. Stability increases from blue to red: values of ^cr on 
lines of equal stability are indicated in the graph. While in both cases stability is optimized for weak self-activation r, they differ 
in the spatial position of the localized front as measured by the value of Mo- (c) For sigmoidal profiles, small Hill coefficients 
n for self-activation are optimal for front stability. Parameters for plots (a-c) were ^ = 100, and k — 0.2. (d) To study if 
regulation of an exponential signal is biologically beneficial we determined the optimal stabilities which can be achieved for a 
front localized at a specific position. For each go there are parameters r, k and Mq such that the linear stability a is maximal. 
Parameters were n = 5, 2 < r < 6, 0.1 < A; < 1. The plot shows the corresponding optimal values for a for exponential (dashed 
line) and sigmoidal external profiles (solid line, m indicated in the graph). Sigmoidal gradients are generally more stable and 
in addition allow stable localization of fronts in a significant distance from the gradients source at a; = 0. 
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Under which conditions is wave localization possible 
and robust? In Drosophila, the parameters r, uq and n 
are of special importance as they are main determinants 
of the gene regulation network |[13| . The wave localizes if 
there is a bistable region in the bifurcation diagram, i.e. 
if, for some x, the reaction term in Eq. ([ij has three real 
roots. Such values of x exist if the maximum value of the 
derivative of Rl^" (u) — m is greater than zero. We obtain 
as an approximate expression for the phase boundaries 
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and ^ ^ 5 + r ■ Figure [ijd) shows that the range of 
allowed parameters grows with n. For large n, the phase 
boundaries are well approximated hy^<^<^ + K 
This translates to gqX kr, i.e. for front localization the 
overall degradation rate at the threshold must be of the 
same order as the maximum production rate. 

To be stable against extrinsic perturbations the front 
should both relax back quickly into its equilibrium po- 
sition and be insensitive to perturbations in the driving 
signal M{x). Since a high relaxation rate implies that a 
front can follow changes in the signal quickly, the two sta- 
bility criteria seem to be somewhat at odds. However, as 
shown below, they are in full accordance with the latter 
being less restrictive. 

The relaxation rate of the front back into its equilib- 
rium position go can be assessed within the framework of 
a linear stability analysis. Mathematically, this is given 
by expanding Eq. (Isl atqo: c(q) = -a{q-qo)+0{q-qo)'^ , 



where a = — dqc{q)\^^^^. The quantity a measures the 
stability of the fixed point qg, such that large values of a 
correspond to a stably localized front. We find 



dM(g)AV{M{q)) ■ dgMiq) 



jrjd.Uix-qWdx 



(5) 



9=90 



revealing that extrinsic stability is determined by three 
factors: In the numerator, the first factor describes how 
sensitively the potential difference of the stable states 
depends on the external source. The second factor, /i = 
I dM{q)/dq 1^^, gives the steepness of the external profile 
at the localization position. While, therefore, a steeper 
source profile enhances front stability, the steepness of 
the front profile (denominator) has the opposite effect. 
The reason simply is that steeper fronts move slower and 
therefore also relax back more slowly, cf. Eq. ([3|. 

Figure |2] shows the results of the numerical evaluation 
of a for both types of external sources; analytical re- 
sults are given in the Supplemental Material. For both 
types of gradients we find that the localized wave front 
is most stable if r is small, i. e if self-activation is weak 
or birth rates are low compared to the strength of the 
external source [Figs. [2] (a) and (b)]. This can mainly 
be attributed to a decreased front steepness: reducing 
self-activation relative to external activation decreases 
the distance between the fixed points and thereby 
the steepness of the wave front. The front's stability is 
further optimized if it is localized at the steepest po- 
sition of the external signal. For signals with a sig- 
moidal profile, this corresponds to Mq « 1/2, and with 
Mq ss Uq — r/2 in dimensionless form, it implies a re- 
lation between the degradation rate and the activation 
rates: aoX = {kr + ku)/'^- Similarly, for an exponential 
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profile with Mq = 1 one finds gqX = kr/2 + Um- 

How does cooperative binding influence stability? 
Since cooperativity in the kinetics of the external source 
implies a steeper sigmoidal profile, large values for the 
Hill coefficient m increase the front's stability; see also 
the explicit expression for a in the Supplemental Mate- 
rial. Conversely, we find that stability is optimized for 
small values of n, i e a low degree of cooperativity in the 
self-activation reaction [Fig. |2] (c)] j^T]. This somewhat 
counterintuitive result can be attributed to a less steep 
front profile for small n; see Supplemental Material. Ex- 
perimental data for Hunchback indeed indicate that the 
Hill coefficient n for self-activation is rather low |13| E5] . 
Figure |2] (d) shows that stability for sigmoidal external 
gradients is, all other things being equal, generally higher 
than for exponential gradients. This implies that regu- 
lating an external positional signal is advantageous to the 
front's stability, since in this case the non- linear ampli- 
fication of the signal makes it possible to create a steep 
signal even far away from the origin. 

To ensure stable localization, the front must also be 
robust against perturbations in M{x). Specifically, it's 
position (7o should only weakly depend on the local signal 
strength: \dq[){M)/dM\j^^j^ <C 1. This condition is equiv- 
alent to a steep source profile, ji — \dM{q) / dq\^^^ ^ 1, 
and hence in full accordance with a large relaxation rate 
a. It is, however, less restrictive since it is indifferent 
to changes in parameters that mainly affect the shape of 
the front, e.g. the rate of self-activation r and the Hill 
coefficient n; see Supplemental Material. 

In many applications the front serves as a signal for 
further downstream processes, e.g. to determine stripe- 
like patterning of the Drosophila embryo |29l I5U] . In 
those instances it is also important that a front is not 
only stable against perturbations, but also sharply dis- 
tinguishes between active and inactive regions. This re- 
quires a steep front which is generally obtained if self- 
activation is strong compared to external activation and, 
to a lesser degree, if binding cooperativity is strong; see 
Supplemental Material. Sharp fronts, however, are sus- 
ceptible to extrinsic ffuctuations and one has to sacrifice 
front stability for the precision of the transmitted signal. 

Intrinsic noise resulting from small copy number fiuc- 
tuations also affects the stability of the localized wave 
front. In this case, stability can be measured in terms of 
the ratio D/Df between the individual particle's and the 
front's diffusion constants. The latter can be calculated 
following the steps outlined in Ref. |17| : 
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C^dx[^{U'n{U) + U{U'' 



(6) 



where h{U) = R^^ {U) + M (x) + U, and U denotes the 
stationary solution. Generally, the front's diffusion con- 
stant is smaller than the particle's diffusion constant by 
a factor N, which corresponds to the typical number of 



particles in the front region. The integral in the numera- 
tor gives the maximum steepness of the front. Hence, as 
opposed to extrinsic stability, intrinsic stability is opti- 
mal for steep fronts. Shallow fronts are prone to stochas- 
tic switching, as the entropy barrier between the stable 
states is reduced in the front region. The terms in the de- 
nominator account for the reaction and diffusion noise. 
In contrast to extrinsic stability, we here find that the 
front is most robust against fiuctuations for strong self- 
activation r. The reason for this is that, as r determines 
the amount of reactions necessary to locally switch be- 
tween the stable states, the rate of stochastic switching 
decreases for large r. Explicit analytical results can be 
found in the Supplemental Material. 

In conclusion, we identified conditions optimizing the 
stability and robustness of localized wave fronts for differ- 
ent types of perturbations. We find that increasing coop- 
erativity in self-activation broadens the parameter regime 
where wave localization becomes possible and thereby in- 
creases the robustness of the localization mechanism. In- 
terestingly, there is a tradeoff between the stability of 
the wave front to extrinsic and intrinsic perturbations. 
While weak self-activation or low birth rates enhance the 
stability with respect to extrinsic perturbations, stochas- 
tic defocussing is minimized for strong self-activation. 
The latter also increases the spatial precision of the sig- 
nal transmitted by the front to downstream processes. 
Moreover, we showed that processing input from exter- 
nal sources with a cooperative gene activation mechanism 
generally enhances the front's stability even far away 
from the source. Surprisingly, while cooperativity in ex- 
ternal activation increases the front's stability with re- 
spect to extrinsic perturbations, the opposite holds true 
for self-activation. 

The conffict between intrinsic and extrinsic stability 
affects, for example, the design of gene circuits in de- 
velopmental systems. Our results suggest different de- 
sign principles depending on the particle number. If the 
number of involved particles is large, intrinsic noise is 
irrelevant. Then the parameters of the genetic network 
may be optimized for robustness against external per- 
turbations which is achieved by weak self-activation and 
strong cooperativity in external activation. Conversely, 
if particle numbers are low, robustness against intrin- 
sic noise requires strong and cooperative self-activation. 
To also safeguard against external perturbation then re- 
quires additional mechanisms beyond those included in 
our simplified model. We expect these general results to 
be important guiding principles in the context of biolog- 
ical pattern forming systems, such as cell polarization or 
the segmentation of embryos. 
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